A local density functional for the short-range part of the electron-electron interaction 
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Motivated by recent suggestions -to split the electron-electron interaction into a short-range part, 
to be treated within the density functional theory, and a long-range part, to be handled by other 
techniques- we compute, with a diffusion Monte Carlo method, the ground-state energy of a uniform 
electron gas with a modified, short-range-only electron-electron interaction erfc(/ir)/r, for different 
values of the cutoff parameter /i and of the electron density. After deriving some exact limits, we 
propose an analytic representation of the correlation energy which accurately fits our Monte Carlo 
data and also includes, by construction, these exact limits, thus providing a reliable "short-range 
local-density functional" . 
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INTRODUCTION AND SUMMARY OF 
RESULTS 



Density functional theorjji*^ (DFT) is nowadays the 
most widely used method for electronic structure cal- 
culations, in both condensed matter physics and quan- 
tum chemistry, thanks to the combination of low com- 
putational cost and remarkable accuracy for a wide va- 
riety of chemical bonds and solid-state systems. There 
are, however, notable exceptions to such an accuracy. 
For example, even the best available approximations of 
the exchange-correlation functional, the key ingredient 
of the DFT, fail to recover long-range van der Waals 
interactions^*^ are not completely safe for the descrip- 
tion of the hydrogen bondi and have intrinsic problems 
with situations of near degeneracy (when two sets of or- 
bitals happen to have very close energies) i^ More gen- 
erally, the "chemical accuracy" (the accuracy needed to 
predict the rates of chemical reactions) has not yet been 
reached. For all these reasons the search for new approxi- 
mate functionals, or even new ways of exploiting the basic 
ideas and advantages of the DFT, is very activeiSiii^iSiLS 

In this context several author s^iSiiSiii have suggested 
to split the electron-electron interaction into a short- 
range part, to be treated within the DFT, and a long- 
range part, to be handled by other techniques. The 
motivation behind these "mixed schemes" is that the 
DFT, even in the simplest local-density approximation 
(LDA), provides an accurate description of the short- 
range electron-electron repulsion,-^ while other tech- 
niques which give a poor description of short-range prop- 
erties, like the configuration interaction (CI) method or 
the random-phase approximation (RPA)fiii^ can, in- 
stead, accurately capture long-range correlation effects. 

Of course there is no unique way to split the Coulomb 
potential 1/r into a short-range (SR) and a long-range 




FIG. 1: Splitting of the Coulomb electron-electron repulsion 
Wee = 1/r into a short-range (SR) part and a long-range (LR) 
part, according to Eqs. Q-lO, with ^ — 1. 
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have been already used for this purpos o^'?'^^ (see Fig.^, 
and we stick to this choice, which yields analytic matrix 
elements for both gaussians and plane waves, i.e., the 
most common basis functions in quantum chemistry and 
solid-state physics, respectively. This form still leaves 
room for some arbitrariness: the choice of the most con- 
venient cutoff parameter fi, which may be different for 
different "mixed schemes" . 

The combination of a short-range DFT calculation 
and a different treatment of the long-range part of the 
electron-electron interaction can be founded on a rigorous 
basis through the adiabatic connection formalism. ^'^'^^^ 
Depending on the specific problem addressed (van der 
Waals forces, near-degeneracy,...), and thus on the par- 
ticular approach to the long-range part of the electron- 
electron interaction, different "mixed schemes" have been 
proposedi2i2ii£ But in all of them, as in standard DFT, 



a crucial role is played by the exchange-correlation func- 
tional, which now must be built for a modi fied el ectron- 
electron interaction. The schemes of Refs. llCJllI need a 
pure short-range functional, E^f-[n\, whose LDA version 
is given by 



e: 



SR, LDAr 



= / n{r)exc{n{r),fi)dr, 



(3) 



where exc{n, fj.) is the exchange-correlation energy per 
electron of a uniform gas of density n interacting with 
a short-range potential like Eq. QJ. The value of n in 
Eq. Q can be either a constant, or, possibly, a conve- 
nient function of the density, /i = /i(n(r))ii& The local 
functional excin, fi) is the quantity which we provide in 
this paper. We start from a jellium-like hamiltonian (in 
Hartree atomic units used throughout this work) 

H = -lf:K+V!I' + KT + V,l^, (4) 



where Vj.^ is the modified electron-electron interaction 



N 



ee r) / , 



1 Y^ erfc(^|rj -Tjl) 



(5) 



Vj,^ is, accordingly, the interaction between the electrons 
and a rigid, positive, uniform background of density n = 

(4^r3/3)-i 



N 



KT-~nJ2 d^- 



aSR 



rfc(^|ri -x|) 



(6) 



and Vf,!^ is the corresponding background-background in- 



teraction 



T^r-yMx dx' , 



, erfc(/^|x — x'l) 



(7) 



First we calculate the ground-state energy per electron of 
this model hamiltonian, as a function of the density pa- 
rameter Ts and of the parameter /i, with a diffusion Monte 
Carlo method (Sec. ^J. Then we derive the asymptotic 
behaviors of the correlation energy €c{rs,fJ,) (Sec. IIII|I . 
On these grounds we finally (Sec. IIV|I present a conve- 
nient analytic parametrization of the correlation energy, 
thus following in the footsteps from quantum simulations 

of the regular jellium model to the best available LDA 
functionalsii^ii8,i9,20 



II. DMC CALCULATION OF THE 
GROUND-STATE ENERGY 

A local density functional for the short-range potential 
of Eqs. ©-iQ should recover the Ceperley-Alderii (CA) 
correlation energy for /i — > 0. In this Section we outline 
the implications of this condition on the technical aspects 



of our calculation, which is in all respects a standard 
application of the diffusion Monte Carlo method in the 
Fixed Node approximation (FN-DMC)^ 

The FN-DMC method gives the energy EpN of the 
lowest-lying Fermionic eigenstate of the Hamiltonian 
which has the same nodes as the chosen trial function 
^T- The error in EpN is variational, and it vanishes 
as the nodal structure of '^t approaches the (unknown) 
nodal structure of the exact ground state. The simplest 
choice for the trial function of a homogeneous fluid^^ is 
the Jastrow-Slater form, '^t{R) = J{R)D{R), where the 
symmetric Jastrow factor J{R) = exp[— J2i<j ^('^y )] d^" 
scribes pair correlations, and D is the product of one 
Slater determinant of plane waves (PW) for each spin 
component {R denotes the coordinates of all the parti- 
cles). A better nodal structure is provided by the so- 
called backflow (BF) wave function. •^^ 

The method used in Ref. J/7 is in principle exact: it 
starts from the FN solution and then it performs a "nodal 
relaxation" , whereby the energy converges to the ex- 
act ground-state result. However, this second process 
is accompanied by an increasing statistical noise, which 
may hinder full convergence of the results. In practice, 
the results of Ref. [ij are between the FN energies re- 
cently calculated with PW and BF nodesSS, and actually 
somewhat closer to the former. Since, on one hand, BF 
calculations are considerably more demanding, and, on 
the other, the most widely used local-density functionals 
are constructed to fit the quantum Monte Carlo results 
of Ref. IttI we choose to stick to the simple trial func- 
tion with Slater determinants of plane waves. In this 
way our "short-range local-density functional" will con- 
tinuously merge into the Ceperley-Alderiii-based local- 
density functionals as ^ ^ 0. 

All the other errors in the simulation can be controlled 
and eliminated. It is easy to ensure that the biases due 
to a finite time step and a finite population of walkers- 
are much smaller than the statistical uncertainty of the 
CA results, which we set as our target precision. The 
number extrapolation is more delicate. We simulate N 
particles in a cubic box with periodic boundary condi- 
tions, interacting via the potential of Eq. |^. Since for 
very small values of fj, we rely on the analytic asymp- 
totic behavior described in Sec. IIIII the only simulations 
we need to do will deal with really short-range poten- 
tials, which we may safely treat using the minimum image 
conventioniS^ The dependence of the energy on the num- 
ber of particles is determined with the Variational Monte 
Carlo (VMC) method, which calculates the expectation 
value of the Hamiltonian operator on the trial wave func- 
tion and is cheaper than DMC. For several values of N 
(namely 38, 54, 66, 114, 162), we use VMC to calculate 
(i) the variational energy Ey (after optimization of the 
Jastrow factor), and (ii) the Hartree-Fock energy Ehf, 
which corresponds to J = 1. For each value of r^ and 
/i, the resulting estimate of the correlation energy per 
electron, ec = {Ey — Ehf)/N, is fitted to the following 
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FIG. 2: Correlation energy per particle of the short-range 
interacting gas ai fi — 0.5 rs = 1, for different numbers of 
particles A'^. The fitting function of Eq. ||HJ (line) favorably 
compares with the VMC data (dots). 



form: 

ecirs,n;N) == edr s, ^J.; oo) + a{rs, fJ,)[T (oo) - T{N)] 
+b{rs,^i)/N. (8) 

Here T{N) is the kinetic energy of N non-interacting 
electrons at rs — 1, and a(rs,fi), b{rs,ij,) and the cor- 
relation energy in the thermodynamic limit, edrs, fJ.',oo) 
are fitting parameters. The size dependence of the VMC 
result for the correlation energy is shown in Fig.|21for the 
case where it is largest (small r^ and small fj,). We point 
out that the simple functional guess of Eq. ||SJ) (solid line) 
accurately models the size dependence of the VMC data 
which, although on a small energy scale, are still far from 
a smooth dependence (dots with error bars). Our final 
result for the correlation energy is obtained by adding 
the infinite-size extrapolation obtained from Eq. IjHJl to 
the result of a single DMC simulation with N = 54. 



III. ASYMPTOTIC BEHAVIORS 



In this section we derive some limiting behaviors of 
the correlation energy edrs, /i), which will be used for its 
parametrization in Sec. IIVI The detailed study carried 
out here can be also of interest for the choice of a density- 
dependent n parameter in the mixed schemes of Refs. 9_ 
andE 

We consider two different regimes: when our system 
approaches the standard jellium model (i.e., full inter- 
action l/r), and when it approaches the noninteracting 
Fermi gas. In the first case (Secs. llTT^ and llIIB|l we find 
that the correlation energy is a function of the scaled 
variable /iy^, while in the second case (Sec. IIII Cll the 
relevant scaled variable is fir^. 



Since for small /i 
erfc(/ia;) 1 

X X 



2^ 



2 x^ 



:M' + 0(A*'), 



(9) 



if we fix the density and let the parameter /x approach 
zero, we can write 



H = H, 



where 



Coul 

N 
2 /I 



^iH^'^+^l''H^^^ + o{^i^), 



(10) 



(11) 



-\ dx dx' |x- X 



^E 



dx Ix — r,- 



(12) 



In Eq. H1U|) . and in the rest of this paper, the suffix "Coul" 
indicates quantities of the standard uniform electron gas 
(jellium), with Coulomb interaction l/r. Thus, for small 
/i we are perturbing the jellium model. 



Ei^,) -- 


= Ecoul + ^iE^^^+^i^E^^^ + ... 


(13) 


*(m) = 


= *Cou1 + M*^''+m'*^'^ + - 


(14) 



Since H^^^^ is a constant, we immediately find E^^^^ = 
i/(i) = N/y^ and *(^) = 0, which, combined with 
i/(2) = 0, also gives ^f^) ^ q ^nd S'^) ^ q ^hen E'-^^ is 
simply 



ij(3) 



(* 



Coul I 



H(3)|^ 



Coul/ 



(15) 



and can be easily evaluated, since it is related to the 
plasma oscillationf24 



£;(3) 



N 



-N- 



3/2 



(16) 



Eqs. (|13|l - 115|l hold because the expectation values of 
H^^^ and H^^^ on ^coui exist, as it will be more explicitly 
shown in Eqs. l(THl) - (|n|l . 

Taking the energy per particle e = E/N, and dividing 
it into the non- interacting kinetic part tg — j^kp and the 
exchange-correlation contribution e^c, we then have the 
small- fj, expansion 



ixc{rs,iJ- -> 0) 



a/tt 



,3/2 

/3^ 



M^ + 0(a.4). (17) 



The same result can be obtained by differentiation of 
E{ii) with respect to /i and by using the Helmann- 
Feynmann theorem, which leads to the exact expression 
(see also Ref.ii): 



djj, 



dss2e-''''-'^'[<7,,(s,r,,^)-l], (18) 



where s ~ r/rg, and gxcis,rs, fJ.) is the pair-distribution 
function^^'^^:— corresponding to the Hamiltonian of 
Eq. Q). The evaluation of Eq. ifT^ at /i ~ 0, imme- 
diately gives the first-order result, l/v^- Higher-order 
derivatives of exc at /x = can be obtained by further 
differentiating Eq. H18|l . provided that the conditions for 
differentiation under the integral sign are fulfilled. Since 
^(1) ^ ^(2) ^ implies dgxc{s,rs,fi)/d^i\f,^o = and 
d"^ gxc{s , rg , fi) / dfi'^l^^o — 0, the possibility of extracting 
the second and third derivatives of e^c at /i = from 
Eq. l(TF|) depends on whether the integrals 



take the limit /i ^- of Eqs. (|24() and H25() we recover 
the first-order result in Eqs. (|22|l and H23|l . However, 
higher-order derivatives at /i = of e^ and ec cannot be 
obtained by differentiating Eqs. (|24|) and (|25|) . This is 
due to the long-range tail of gx{s) — 1 and g^°'^^{s,rs): 
when s ^ oo they both approach zero as^LSLSS l/s^. 
Thus, integrals of the kind J^ s'^[gx ~ 1] and J^ s^g^°'^^ 
diverge. The long-range tails of gx{s) — 1 and g^°^^{s, Tg) 
exactly cancel2Si2L2a in g^°'^\s,rs) - 1. This is why, at 
small /i, both ex and Ec have terms ex /i^ which cancel 
out in exc- 



dss^[g:^ri-s,rs)-l] 



(19) 



with n 



4 and n 



This is the case, 



6 exist 
since gx. 

oscillation-averaged parliSS goes to zero as^^j^l 1/s* when 
s 



Couij-^^^^-j — 1 is a well-behaved function whose 



oo. We thus obtain from Eq. ((TF|) 



d^exc 



dfi^ 



Ai=0 



Ai=0 







(20) 



6 



dss'[g^r\s,r,)~l] 



r.3/2 



r gU,p 



r:-, (21) 



in agreement with Eq. (|17|l . We see that since 5^°"'(s -^ 
cxD, r^) — 1 (X 1/s® no further information can be extracted 
from Eq. p8|) . or, equivalently, by going further with the 
expansion of Eq. Hl()|l . 

One can divide exc into its exchange and correlation 
parts, exc = ^x + Cc- The exchange energy ex has been 
calculated by Savin(Si and is reported in Appendix^ Its 
small-/i expansion is 



exirs,iJ.^ 0) 



,Coul 



{rs) + 



M 



3ar 



2tt 



s 2 



-0{^l% (22) 



where a — (4/97r)^/^. The /i ~+ behavior of ec = 
Exc - £x, is then 



ecirs,iJ.^ 0) 



e"V.) + ^"''^ ■■' 



277 



-M 



3/2 



M' + 0(/i^)- 

(23) 

Notice that if we divide the pair-distribution function gxc 
into its exchange and correlation parts, gxc{s,rs, ^j) — 
9x{s) + gcis,rs,fj.), we have 



dex 

dec 
9/i 



^/ dss2e-^^'-=^^[g.(.)-l], (24) 
^ Jo 

dss^e-^"'-"'"gc{s,rs,fi). (25) 



^71" Jo 



(This follows directly from the Hellmann-Feynmann the- 
orem and from the fact that gx corresponds to the non- 
interacting gas and thus does not depend on /i.) If we 



B. Finite /i, and r^ ^ 

If we use the relevant scaled units s — r/r^ and we let 
Ts approach zero, the potential has the expansion 



1 ericiprgs) 



1 



2fi 



-— r 

3 ^/^ ^ 



which has the Coulomb interaction as leading term. We 
are thus approaching again the jellium model, so that 
Eq. H23() is also valid for finite /i and r^ -^ 0. 

In Eq. H23|) the relevant scaled variable is fJ.^/r^- This 
can be understood in the following way. The Coulomb 
gas presents screening effects at lenghts r > l/gxF oc 
Y^, where gxF is the Thomas-Fermi wave vector. Since 
the erfc function amounts to some sort of artificial screen- 



ing at lenghts 



> 



Ai 



the Thomas-Fermi screening ap- 



pears, exactly as in the Coulomb gas, when y^ ^ /i 



urs » 1 



When /i -^ oo, the potential terms of Eqs. ®-0 



rapidly vanish (V 



SR 



-fi r' 



). In this regime we can 



treat the whole potential as a perturbation to the non- 
interacting Fermi gas. 

The first-order (in the potential) correction to the non- 
interacting energy tg gives ex of Appendix^ The second- 
order term can be computed by standard Rayleigh- 
Schrodinger perturbation theory 



N 



-y- 

nj^O 



V^^\n){n\V''"-\0) 

En — Eq 



(26) 



As in the case of jellium, e^^^ is the sum of a direct term 
and of a second-order exchange ternij^ which in Fourier 
space read 
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rfq 



dq 

1 



dk 



|k+q|>l 



dp ■ 



dk / dp 

|k+q|>l -/ip+q^l 



|p+q|>i 1 
1-e- 



[l-k)0{l-p) 
q-(k + p) 



r 



|q + k + p|' 



d{l-k)9{l-p) 

^2 + q . (k + p) • 



(27) 
(28) 



Here all the momenta are expressed in units of kp — 
{ars)~^, and 9{x) is the Heaviside step function. Now, 
consider the case fj^rg — s- cx) and divide the integral over 
q in Eqs. (|27(l and H28(l into two parts: 



dq — dq 

Jo 



dq, with 1 ^ qi ^ /irs 



/O JO Jqi 

In the first part, when q G [0, gi], we can write 



1 _ g-g fef/V 



1 _ g-lq+k+pl^^/V 



4/^2 ' 

|q + k- 



p|2fc 



21,2 



4^2 



(since gi ^ /zr^, and the integrals of Eqs. (|?7jl and l(^ 
are restricted to the domain |k| < 1, |p| < 1). Equa- 
tions H27|l and (|28|) then reduce to integrals of the same 
kind, which can be summed to yield 



-|q|<9i 



-3 



327r5 \An 



kp 



dq / dk 

|q|<gi ^|k+q|>l 

l~k)0{l-p) 



/|p+qi>i r + q-(k + p) 



i.e., they give a term which vanishes as (/ir^)^^. In the 
second part, q G [qi,oo), having chosen gi ^ 1, we can 
write 



|q + k + p|' 



<z^ 



dk / dp^^L^W^ 

ik+qi>i -'ip+qi>i r + q-(k + p) 



Equations (|27l) and H28() again reduce to integrals of the 
same kind, which can be summed to yield 




.(2) 

1q|>gi 37^2 



(29) 



The right-hand side of Eq. H29I) can be evaluated analyti- 
cally and then expanded for iiVg —> oo. Its leading term is 
(correctly) independent of gi and equals — /'^ (/^'"s)"'^- 
We thus have 
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with A « #^ 
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FIG. 3: Correlation energy of the short-range interacting gas 
as a function of /ir-s, for different densities. Our fitting func- 
tion (lines) is compared with our DMC data (dots). The error 
bars are comparable with the symbol sizes. 



Since the perturbation series expansion whose second- 
order term corresponds to Eq. (|26() is done with respect 
to the whole potential V^^ and not with respect to the 
parameter /i, higher-order terms could also contribute to 
the value of A. For this reason, in our parametrization 
of ec(rs, C) ^ is left as a free parameter, to be optimized 
with a fit on the DMC data. We expect to find a value of 
A of the same order of the one estimated with Eq. H26|l . 
since the potential V^^ vanishes very rapidly as /i — > oo, 
so that the higher-order-term contribution to A should 
be small. 



IV. ANALYTIC REPRESENTATION OF THE 
CORRELATION ENERGY 



An accurate and simple analytic representation of the 
correlation energy ec{rs,p) can be obtained by a Fade 
form which interpolates between the limiting behaviors 
given by our Eqs. (|23|l and (|30|l . and contains some free 
parameters to fit our DMC data. We write 



:(?'s,M) 



,Coul 



(r,)[l-h6i(r,)/i] 



1 + fei(r,)/x -f- 62(r,)/x2 -f- 63(r,)^3 + hi{r,)^i^ ' 

(31) 



where 



61 = 

62 - 
64 = 



63- 
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3/2 , 

Ts' 1 


V3^ e?°"i 
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. ^Coul ' s 

-he, -J, 



(32) 
(33) 
(34) 



and e^°"'(rs) is one of the standard 
parametrizations^^'-'^^'^'i of the correlation energy of 
the unpolarized jeUium. Here we used the paranietriza- 
tion of Perdew and Wang.^*' The two parameters &3(?'s) 
and A are fixed by a two-dimensional (rj , fi) best fit to 
our DMC data. We find: 



A = 0.03579. 



(35) 
(36) 



This fit yields a reduced x^ of 2.7. In Fig. 01 we show our 
DMC data together with the fitting function for different 
values of rg. Notice that our analytic €c{rs,^i) does not 
break down at high (rs — *■ 0) or low (r^ -^ 00) densities, 
being constrained by exact behaviors. 



will be the generalization of our functional to the spin- 
polarized case. 

What we obtained in this paper is not the only pos- 
sible short-range local-density functional of interest to 
"mixed schemes" . In some of them.^ the DFT treatment 
of the short-range part is handled through another func- 
tional E^, [n] , defined as the difference between the stan- 
dard exchange-correlation energy functional (correspond- 
ing to the Coulomb interaction) and a long-range-only 
functional 



-SR 



E,M^E,c[n]-E^'^\n]. 



(37) 



Another direction of future work will thus be the study of 
the uniform electron gas with a long-range-only interac- 
tion of the form of Eq. (0) , and, possibly, other modified 
interactions proposed in the same spiriti^ 
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V. CONCLUSIONS AND PERSPECTIVES 

We have presented a comprehensive numerical and an- 
alytic study of the ground-state energy of a (spin unpo- 
larized) uniform electron gas with modified, short-range- 
only electron-electron interaction erfc(/ir)/r, as a func- 
tion of the cutoff parameter /i and of the electronic den- 
sity. Our chief goal has been the publication, in a con- 
venient form for application, of a reliable local density 
functional for the correlation energy of this model system, 
which (i) fits the results of our quantum Monte Carlo sim- 
ulations and (ii) automatically incorporates exact limits. 
Such a functional is a crucial ingredient for some recently 
proposed "mixed schemes" , which exploit the DFT only 
for the short-range part of the electron-electron interac- 
tion. In this context the natural extension of this study 



APPENDIX A: EXCHANGE ENERGY 

The exchange energy corresponding to the Hamilto- 
nian Q has been calculated by Savin in Ref. 1^ and is 
equal to 



ej;(rs,/i) = kpl - - a 

n a 



y^cii ( — ) - 3a + 4a^ 



+ (2a-4a=^)e-i/*"' 



(Al) 



with a = /i/(2fci?). The exchange energy thus satisfies 
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